###与其他不同的是只利用了2018sci附件4的TWAS数据，而且没有用TWAS不显著的数据
setwd("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/")
library(ChIPpeakAnno)
library(meta)

col_names=c("TWAS.db","con.file.name","case_in_region","case_not","con_in_region","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=8))
names(result)=col_names
fn=c("SCZ.TWAS","ASD.TWAS","BD.SCZ")
for(i in 1:length(fn)){
file=read.xlsx("aat8127_Table_S4.xlsx",colNames = T,sheet = fn[i])
file$TWAS.P=as.numeric(file$TWAS.P)
file_sig=file[file$TWAS.P<0.05,]
file_sig=file_sig[!is.na(file_sig$CHR),]
file_sig$unitID=paste0("chr",file_sig$CHR,":",file_sig$start,":",file_sig$stop)
file_sig=file_sig[!duplicated(file_sig$unitID),]

file_t=data.frame(chr=paste0("chr",file_sig$CHR),start=format(as.numeric(file_sig$start),scientific = F),end=format(as.numeric(file_sig$stop),scientific = F))
bed_file=toGRanges(file_t,format="BED",header=TRUE)

ol1=findOverlapsOfPeaks(bed_case,bed_file)
tcase=as.data.frame(ol1$peaklist$`bed_case///bed_file`)
for(j in 1:4){
con.file=c("nobias_AShM.txt","nobias_AShM_total_reads10.txt","nobias_AShM_p.5.txt","nobias_AShM_total_reads10_p.5.txt")
con=read.table(con.file[j],head=F,sep="\t")
con_t=data.frame(chr=con[,1],start=as.numeric(con[,2]),end=as.numeric(con[,3])+1)
bed_con=toGRanges(con_t,format="BED",header=TRUE)

ol2=findOverlapsOfPeaks(bed_con,bed_file)
tcon=as.data.frame(ol2$peaklist$`bed_con///bed_file`)

result_tmp=data.frame(matrix(NA,1,ncol=8))
names(result_tmp)=col_names

a=dim(tcase)[1]
b=727-a
c=dim(tcon)[1]
d=dim(con_t)[1]-c
result_tmp[,1]=fn[i]
result_tmp[,2]=con.file[j]
result_tmp[,3]=a
result_tmp[,4]=b
result_tmp[,5]=c
result_tmp[,6]=d
result_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
result_tmp[,8]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
result=rbind(result,result_tmp)
}
}
metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=result,sm="OR",studlab = TWAS.db) 
forest(metaor3)